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ABSTRACT 

The  occurrence  of  independent  events  at  random  in  the 
plane,  i.e.  the  formation  of  a  planar  point  process,  is 
discussed.   Both  homogeneous  and  nonhomogeneous  processes 
are  considered.   A  specific  functional  form  for  the  parameter 
in  a  nonhomogeneous  planar  Poisson  process  is  used  to 
illustrate  the  development  of  test  and  parameter  estimation 
techniques.   The  problem  finds  application  in  the  description 
of  biological  phenomena  as  well  as  in  search  and  detection 
problems . 
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I.   INTRODUCTION 

Many  problems  arising  naturally  in  a  physical  sense  are 
often  so  complex  that  the  identification  and  description  of 
underlying  mechanisms  must  use  the  tools  of  probability  and 
statistics.  Some  of  the  reasons  leading  to  the  requirement 
of  using  these  tools  are: 

(i)    the  data  base  may  be  so  large  or  complex  as  to 

preclude  identification  of  any  driving  mechanism 
without  recourse  to  statistical  analysis; 

(ii)   if  identifiable,  the  mechanisms  may  be  inherently 
probabilistic;  or 

(iii)  if  identifiable  and  deterministic ,  the  governing 
law  which  the  mechanisms  obey  may  be  unknown. 

This  paper  is  concerned  with  the  use  of  statistics  in 
the  identification  and  mathematical  description  of  the  sp?  ial 
distribution  of  events  (occurrences).   Included  is  the  det  c- 
tion  and  estimation  of  parameters  which  influence  the 
description  of  this  distribution. 

The  area  of  concern  here  is  a  departure  from  those  sta- 
tistical methods  which  have  been  developed  to  detect  the 
effect  of  varying  a  controlled  segment  of  the  underlying 
mechanism.   Among  those  methods  would  be  the  design  of  exper- 
iments, regression  analysis,  time  series  analysis,  and 
analysis  of  variance.   One  goal  of  such  analysis  is  to 
hopefully  predict  the  advisibility  of  pursuing  some  course 
of  action. 


In  the  basic  model  of  this  paper,  events  are  considered 
to  occur  v;ith  a  Poisson  distribution  in  the  plane.   This 
"is  the  natural  model  for  the  expression  that  'points  are 
distributed  at  random',"  [Fisher,  1972,  p.  1*11].   The  bi- 
variate  Poisson  process  will  be  defined  and  then  developed 
through  the  use  of  partial  differential-difference  equations, 
a  widely  repeated  procedure  in  the  univariate  case  but 
neglected  in  the  bivariate  case. 

Initially  a  homogeneous  Poisson  process  will  be  assumed 
to  control  the  underlying  mechanisms.   Then  trends  will  be 
introduced  by  defining  the  Poisson  parameter  in  such  a  way 
as  to  make  it  be  spatially  dependent.   This  will  be  the  basis 
for  the  definition  of  the  non-homogeneous  Poisson  process. 
Time  inhomogeneity  will  not  be  considered.   Thus,  the  data 
are  assumed  to  be  taken  concurrently,  i.e.,  the  period  of 
observation  is  short  compared  to  any  period  of  change  of 
the  parameters. 

Tests  will  be  developed  to  distinguish  between  homogeneity 
and  non-homogeneity  and  the  method  of  maximum  likelihood 
will  be  used  to  develop  estimates  of  the  parameter  in  the 
homogeneous  case  and  parameters  in  the  non-homogeneous  case. 
In  the  latter  case,  conditional  likelihood  techniques  will 
be  utilized  to  develop  tests  and  estimates.   Throughout, 
testing  and  estimation  procedures  will  be  based  on  a  single 
realization  of  the  process  which  consists  of  the  number  of 
events  observed  and  their  spatial  locations. 


The  problem  of  concern  finds  application  in  the  estima- 
tion of  the  density  of  trees  in  a  forest;  here  one  might  be 
concerned  with  estimating  the  potential  yield  of  lumber  from 
a  given  forest  area  where  inhomogeneities  arise  due  to  soil, 
weather  patterns,  topography  and  other  physical  rea.sons . 

Another  application  might  be  in  naval  search  and  detec- 
tion problems.   For  example,  one  might  be  searching  for  a 
merchant  ship  in  distress  whose  location  is  not  known  exactly 
due  to  failure  of  the  ship's  communication  equipment.   Here 
the  independence  assumptions  of  the  planar  Poisson  process 
may  be  valid,  but  not  the  assumption  of  homogeneity.   In- 
homogeneities of  location  occur  because  of  preferred  sea 
lanes  and  physical  characteristics  of  the  ocean  and 
atmosphere . 


II.   THE  HOMOGENEOUS  POISSON  PROCESS  IN  THE  PLANE  (HPPP) 

A.   GENERAL  DEVELOPMENT 

Consider  a  stochastic  process  of  events  occurring  in 
the  plane  (i.e.,  a  so-called  point  process)  which  is 
characterized  by  the  assumptions 

I.  There  exists  a  finite  positive  constant  X  >   0. 

II.  For  any  integer  k  _>  1  and  any  set  of  non-overlapping 
regions  R.  , '  *  '  ,  R,  with  areas  A,,*"*, A,,  (in  the 
usual  geometric  sense),  the  number  of  events  occurring 
in  any  region  i,  denoted  N(R.),  has  a  Poisson  dis- 
tribution with  parameter  XA.  which  depends  only  on 
the  area  of  the  region,  A.,  and  not  its  shape.   Thus, 

ni 
(XA.)  Xexp(-XA  ) 

prob  {N(R. )  =  n  }  = — — , —  .   (1) 

III.  Further,  N(R..  ),  i  =  l,2,"'*,k,  are  mutually  indepe  - 

dent  in  that  N(R. )  is  not  affected  by  the  occurren  e 

of  events  in  any  other  region  or  in  any  grouping  of 

the  regions,  G,  as  long  as  R.flG  =  0.   Thus 

n.  -XA. 
k   (XA  )  1e    X 

prob(N(R,  )=n,  ,  i  =  l,"*,k}=   n   i— — } (2) 

1    x  1=1      ni* 

Definition  1:   If  a  process  obeys  the  above  assumptions  it  is 
called  a  homogeneous  planar  Poisson  process  (HPPP). 

For  reasons  of  arbitrary  shape  the  above  basic  definitions 
will  suffice.  However,  under  certain  geometrical  assumptions,  ai 


equivalent  definition  for  the  HPPP  can  be  achieved  in  a  man- 
ner similar  to  the  development  of  the  univariate  Poisson 
process  through  the  use  of  partial  differential-difference 
equations.   This  is  useful  for  the  development  of  statisti- 
cal properties  and  will  be  very  important  in  the  developmenl 
of  the  non-homogeneous  process.   Such  a  development  also 
provides  another  phenomenological  approach  to  the  ho.:,  genecus 
Poisson  process,  one  which  might  arise  through  the  struc- 
turing of  a  model  for  instance.   For  illustrative  purposes 
the  following  development  will  be  accomplished  using  rectan- 
gular regions.   Note  that  the  development  is  very  dependent 
on  the  geometry  involved;  hence  developments  with  other 
geometries  (e.g.  circular  regions)  must  proceed  somewhat 
differently. 

The  underlying  assumptions  in  the  differential  equation 
development  will  be 

I'.   There  exists  a  finite  positive  constant  X  >  0. 

II'.  For  any  region  R*  with  incremental  area  AAS  inde- 
pendent of  the  shape  of  the  region  except  possibly 
as  noted  above 

(a)  prob  {no  event  in  R*}  =  1  -  XAA  +  o(AA), 

(b)  prob  {one  event  in  R*}  =  XAA  +  o(AA), 

(c)  prob  {more  than  one  event  in  R*}  =  o(AA), 

t?C  AA ) 
where  "g(AA)  is  o(AA)"  means  lim  fe  ^a   =  0 ,  or 

AA+0    AA 

specifically  in  rectangular  regions  the  limit  as  Ax 

o-(AA) 
or  Ay  or  both  go  to  zero  of  °  .    is  zero. 


III'.   The  occurrence  of  events  In  R*  is  independent 
of  the  occurrence  of  events  in  any  region  R 
where  R*f|R  ■  0- 
It  will  be  shown  that  If.  11 '  and  III1  imply  and  are 
implied  by  I,  II  and  III  so  'chat  the  two  sets  of  assumptions 
are  equivalent  and  hence  the  incremental  assumptions  give 
rise  to  a  HPPP.   Clearly  _  and  I1  are  the  same,  as  are  III 
and  III'.   Also  II  implies  II1  since  by  (1) 


2 
(a)   prob  (N(R*)  =  0}  =  e~XAA  =  1  -  XAA  +  ~   (AA)2 

p 
X    5   2 
■  1  -  XAxAy  +  ~-  Ax~Ay   -  . . . 


-  XAA  +  o(AA), 


with  the  definition  of  o(AA)  given  above.   Also 


(b)   prob  (N(R*)  =  1}  =  XAAe  XAA  =  XAA(1  -  XAA  +  . . . ) 


=  XAA  +  o(AA) 


oo   ,.  .*i  -XAA 

and  (c)   prob  {N(R*)  >  2}  =  I      lAAA;  * =  o(AA). 

1=2       1# 


The  problem  remaining  in  order  to  demonstrate  equivalence 
between  the  two  sets  of  assumptions  is  to  show  that  II' 
implies  II. 

Consider  a  region  R  bounded  by  the  co-ordinate  axes  and 
lines  x  =  X*  and  y  ■  Y*,  with  area  X*Y*.   Now  extend  the 


sides  tc  x  =  X*+Ax  and  y  =  Y*+Ay  (see  Figure  1).   Consider 
the  probability  of  n  events  occurring  in  the  extended 
region,  R»  =  R  +  Pl]L  +  R2  +  R   where: 

(a)  R  has  area  X*Y*, 

(b)  R  has  area  X*Ay, 

(c)  R_  has  area  Y*Ax, 

(d)  R~  has  area  AxAy; 

(a)-(d)  imply  R1  has  area  X*Y*  +  X*Ay  +  Y*Ax  +  AxAy. 
The  assumptions  IT,  II',  and  III1  imply 


prob  {no  event  in  R  }  =  1  -  AX*Ay  +  o(X*Ay), 

prob  {one  event  in  R  }  =  AX*Ay  +  o(X*Ay),  (3) 

prob  {more  than  one  event  in  R.}  =  o(X*Ay); 

prob  {no  event  in  Rp}  =  1  -  AY*Ax  +  o(Y*Ax), 

prob  {one  event  in  Rp}  =  AY*Ax  +  o(Y*Ax)s  (4) 

prob  {more  than  one  event  in  Rp}  =  o(Y*Ax) ; 

and 


prob  {no  event  in  R_}  =  1  -  AAxAy  +  o(AxAy), 

prob  {one  event  in  R_}  =  XAxAy  +  o(AxAy)3  (5) 

prob  {more  than  one  event  in  R_}  =  o(AxAy). 


Moreover,  statements  (3)5  (*Os  and  (5)  are  independent. 

It  Is  noted  that  the  above  equations  may  have  two 
different  interpretations.   For  instance  in  (3),  prob{one 
event  in  RJ  =  X*Ay  +  o(X*Ay)  is  interpreted  to  mean  one  event 
in  a  two-dimensional  process  with  parameter  X  and  area  X*Ay. 


10 
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Figure  1.   The  incremental  increase  of  a 
rectangular  region. 
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However,  another  interpretation  would  be  to  consider  the  one- 
dimensional  (marginal)  process  of  events  projected  onto  the 
y-axis,  in  which  case  the  parameter  is  XX*  and  the  incremental 
interval  has  length  Ay. 

For  notational  convenience,  let  P  (X*  Y*)  denote  the 

j      n   ' 

probability  that  n  events  occur  in  a  region  with  area  X*Y*. 
The  differential-difference  equations  are  written  noting  that 
n  events  may  occur  in  an  extended  region  by  having  n  events 
in  the  unextended  region  and  no  events  in  the  extension, 
n-1  events  in  the  unextended  region  and  one  event  in  the 
extension,  etc.  Hence 

P    (X*+Ax,Y«)    =    P    (X*,Y*)    •    P    (Ax,Y*) 
n  *  n        '  os 

+   Pn    ,(X*,Y*)    •    P,(Ax,Y*) 
n-1  1 

+   Pn_2(X*,Y*)    ■    P2(Ax,Y*)    +    .. . 

=    P    (X*    Y*)[l-XY*Ax]    +    P      , (X*Y*)[XY*Ax]    +    o(Y*Ax). 
n  n-i 

(6) 
Similarly, 

P    (X*,Y*+Ay)    =    P    (X*,Y*)[l-XX*Ay]    +    P      . (X* ,Y* ) [ XX*Ay ]    + 
n  n  n-x 

+   o(X*Ay),  (7) 

and 
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P  (X*+Ax,Y*+Ay)  =  P  (X*,Y*)[l-XY*Ax][l-AX*Ay][l-XAxAy] 

n  n 

+  Pn_1(X*,Y*)[XY*Ax(l-XX*Ay)(l-XAxAy) 

+  XX*Ay(l-XY*Ax)(l-XAxAy) 
+  XAxAy(l-XX*Ay)(l-XY*Ax)] 
+  Pn^2(X*,Y*)[XX*Ay  •  XY*Ax(l-XAxAy) 
+  XX*Ay(l-XY*Ax)XAxAy 
+  (l-XX*Ay)XY*AxXAxAy] 

+  Pn_3(X*,Y*)[X3X*Y*Ax2Ay2]  (8) 

+  o(AxAy)  +  o(X*Ay)  +  o(Y*Ax). 

Interpreting  the  above  equations,  the  third  term  on  the 
right  hand  side  of  (8),  for  example,  states  that  there  can 
be  n  events  in  the  extended  region  R'  if  there  are  n-2 
events  in  R  and  exactly  one  event  in  each  of  any  two  of  t  a 
added  regions.   That  is,  there  can  be  two  events  in  the 
added  regions  R.  ,  R«  and  R_  if  one  occurs  in  each  of  two 
regions  and  none  occurs  in  the  third  region,  i.e.,  one  in 
R,  ,  one  in  Rp  and  none  in  R-.,  etc.   Collecting  all  terms  of 
order  o(AxAy),  o(X*£y)  and  o(Y*Ax),  (8)  reduces  to 
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Pn(X*+Ax,Y*+Ay)=Pn(X*,Y*)[l-AY*Ax-AX*Ay-AAx£y+A2X*Y*£xAy] 
+  Pn_1(X*,Y*)[AY*Ax-2X2X*Y*AxAy+XX*Ay+XAxAy] 

+  Pn_p(X*,Y*)[X2X*Y*AxAy]+o(AxAy)+o(X*Ay)+o(Y*Ax). 

( 8 f ) 

=  Pn(X*,Y*)[l-XY*Ax]+Pn_1(X*,Y*)[XY*Ax]+o(Y*Ax) 
+  Pn(X*,Y*)[l-XX*Ay]+Pn_1(X*,y*)[XX*Ay]+o(X*Ay) 

-  Pn(X*,Y*)-XPn(X*,Y*)AxAy+XPn_1(X*,Y*)AxAy 

+  X2X*Y*[Pn(X*,Y*)-2Pn_1(X*,Y*)+Pn_2(X*,Y^)]AxAy-cr  •;-.•; 

Noting  from  equations  (6)  and  (7)  that  the  firs-  three 

terms  on  the  right-hand  side  of  the  above  equation  ar : 

P  (X*+Ax,Y*)  and  the  next  three  are  P  (X*,Y*  +  Ay)  ar.;. 
n       5  n    3  J 

rewriting  (8')3  the  result  is 

P  (X*+Ax,Y*+Ay)=P  (X*+Ax,Y*)+P  (X*,Y*+Ay)-P  (X*aY»] 
n  n  n    '  n 

-  XPn(X*,Y*)AxAy+XPn_1(X*,Y*)AxAy"  (8") 
+  X2X*Y*[Pn(X*,Y*)-2Pn_1(X*,Y*)+Pn_2(X*,Y*]AxAy 

+  o(AxAy) . 

The  definition  of  the  second  partial  derivative  with 
respect  to  two  variables  is 
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l_{3F(x,y)}=lim  _i_{lim  ?(x+Ax5y+^y)-F(x,y+Ay)  _ 


9y  9X       Ay-O^  Ax-K) 


.  lim  F(x+Ax  yJ^ZCx^  > 
Ax+O       Ax 


Hence,  transposing  the  first  three  terms  of  equation  (8") 
to  the  right  hand  side,  aivLdir>£  V,  AxAy  and  then  taking 
the  double  limit  result"  I:, 


3^P  (X*  Y*)  ? 

n   ___  =  -XP  (X»,Y*  )+XP   -,(X*.Y*)+A^X*Y*[P  (X*  Y*) 

n  n--  n 


8x8y 


-2P   .(X«fY«)+Pn  ,vX*,Y*)]  (9) 


The  solution  to  (9)  (a  p,  r,.  :.*.  differential-difference 
equation)  is  easily  show?  •   be 


P  (X*,Y«)  =  K(XX*Y«)"e>.p(-XX»Y»)  ,   n  =  2,3, (10a) 

n  n ! 


where  K  is  an  arbitrary  constant.   Special  considerations 
are  needed  for  n  =0,1  since  for  these  cases  some  of  the 
terms  in  (8")  and  (9)  are  not  defined.   Rewriting  (8")  and 
(9)  while  concurrently  eliminating  the  proper  terms  leads  to 


P  (X*,Y*)  =   K(XX»Y»)nexp(-XX»Y»)  ,   „.  0.1.2. ... (10b) 


n  *  n 


Since  P  (X*  Y*)  is  a  probability  statement  and  for  any 
n 

given  region  the  number  of  events  in  that  region  must  be 
some  non-negative  integer,  the  constant  K  is  seen  to  be  unity 
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Hence  (10b)  is  equivalent  to  (1)  which  was  to  be  proven. 
Thus  the  two  sets  of  assumptions  imply  the  same  things, 
namely  that  the  number  of  events  in  a  region  has  a  Poisson 
distribution  with  parameter  proportional  to  the  area  of  the 
region  and  independent  of  its  shape  and  the  number  and 
position  of  events  outside  the  region.   Note  that  the 
formulation  excludes  multiple  events,  i.e.,  the  occurrence 
of  two  or  more  events  at  any  point  or  on  any  line  in  a  single 
added  region  such  as  R  in  Figure  1. 

Also  a  similar  derivation  will  go  through  for  circular 
regions  using  polar  co-ordinates,  but  there  are  differences 
in  the  special  properties  of  the  Poisson  process  as  defined 
through  assumptions  I,  II  and  III  in  differently  shaped 
regions.   These  are  discussed  below.   The  differences  in 
the  special  properties  of  the  non-homogenecus  planar  Pols  sen 
processes  as  they  vary  with  different  geometries  are  an 
essential  element  of  the  analysis  of  points  (events)  in 
the  olane. 

B.   TESTING  DATA  FOR  HOMOGENEOUS  PLANAR  POISSON  PROCESS  (.HPPP) 

Given  the  occurrence  and  spatial  location  of  n  events  in 

a  rectangular  region  of  area  X*Y*,  consider  the  problem  of 

determining  whether  or  not  these  points  occur  as  realizations 

of  the  HPPP.   Miles  [1970, p. 89]  has  stated  a  consequence  of 

definition  1  as 

Corollary.   Assume  a  rectangular  region  R.  with  area  A  . 

Given  N(R.)  =  n  and  0  <  A,  <  °°,  the  n  points  are  independently 
1  i 

and  uniformly  distributed  in  R. . 
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Proof:   Let  A.  =  X*Y*  where  R.  is  a  rectangular  region 
bounded  by  the  coordinate  axes,  x  =  X*  and  y  =  Y4: .   Label 
the  n  given  points  in  any  convenient  manner,  i  .g. ,  on  the 
magnitude  of  the  y-component.   Let  (xj),.,  denote  the  i  ] 
labelled  event.   Consider  an  incremental  region  with  area 
dxdy  which  has  the  property:   prob  {exactly  one  event  in  the 

incremental  region  of  area  dxdy}  =  P-,(dx,dy)  =  Xdxdy  +  o(dxdy) 

j. 

Consider  now  n  incremental  recangles  dx.dy.,  i  =  l,"**,n, 
placed  in  R .  .   Ignoring  probabilities  of  o(dx.dy.),  assump- 
tions I,  II  and  III  imply  that  the  joint  probability  that 

x.  i_ 

the  i   event  falls  in  the  incremental  rectangle,  dx.dy., 
i  =  1,  •••,!!,  and  exactly  n  events  occur  altogether  in  X*Y* 
is  given  by 

Xdx1dy1. . .Adxndynexp{-AX*Y*}. 
Restating  in  terms  of  the  density  function, 

f{(x,y)(1),. . .  ,(x,y)(n),n;X}  =  Anexp{-AX*Y*} , 

where  f{...}  is  the  joint  density  of  (x,y),±y    i  -   l,...,n, 

and  the  probability  that  the  number  of  events  in  X*Y*  is  n. 

The  exponential  term  in  the  above  expressions  is  an  approx- 

n 
imation  to  exp(-(AX*Y*  -   I   Adx.dy. )},  i.e.,  represents 

i-1 
the  probability  of  no  events  within  the  region  X*Y*  but 

outside  the  incremental  regions  containing  each  event. 

By  conditioning  on  the  occurrence  of  n  events  in  the 

region  which  are  distributed  Poisson  with  parameter  AX*Y*, 
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«x,y)a),...(x,y).,Jn;X)  =  *"exp(-AX«Y») 

1  '  (    '  (XX*Y*)exp(-XX*Y* 


-,  n  >  1 

) 
n ' 


n!  (11) 


(X*Y*) 


n    ' 


which  is  the  joint  distribution  for  n  bivariate  uniform 
random  variables  ordered  on  one  of  the  random  variables  as 
is  shown  in  Appendix  A.   Note  also  the  independence  of  the 
conditioned  density  from  the  parameter  A,  i.e.  the  random 
variable  N  is  a  sufficient  statistic  for  A. 

As  a  consequence  of  the  above  corollary,  it  is  apparent 
that  if  the  points  of  the  HPPP,  conditioned  on  the  number 
of  events  observed  to  occur,  are  in  fact  ordered  with 
respect  to  the  increasing  magnitude  of  the  y-component,  then 
no  "information"  is  available  about  the  ordering  of  the  x- 
components,  i.e.,  each  of  the  n!  orderings  that  can  be 
induced  on  the  x's  by  the  orderings  on  the  y's  has  probability 
1/n!.   This  is  readily  apparent  since  in  the  bivariate  uniform 
case  the  two  components  were  independently  selected.   Hence, 
if  (x,y),,x'is  determined  by  (x^y/vO*  i.e.  the  points  are 
labelled  by  the  ordered  y-component,  then 


prob(Xk  =  X(j)>  =  ~    J  =  1,2,. ..,n 


where   X,.,.    is   the   j        X     in  magnitude,    and 


prob(X1   =   X,.*,    j  =  l,...,n;    X2=X(k),    k=l , . . . , j-1 , j+1 , . . . n; 

...;    X      =   Xr(M}    =   ^r  (12) 

9      n  (x,)  n! 
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Hence  if  the  x-components  of  the  points  ordered  on  the 
y-components  exhibit  any  natural  ordering  then  the  x-  and 
y-compcnents  have  not  been  independently  selected  and  the 
observed  process  cannot  be  a  HPPP.   This  will  be  the  basis 
for  many  of  the  tests  for  a  HPPP  against  a  non-homogeneous 
planar  Poisson  process  to  be  discussed  later. 

Lemma :   If  the  bivariate  process  is  Poisson  and  the  regions 
are  rectangular,  then  the  projections  of  the  events  onto 
the  coordinate  axes  may  be  shown  to  be  univariate  Poisson. 

Proof:   Consider  the  occurrence  of  events  in  a  rectangular 
region  of  area  X*Y*.   Then  by  III  the  occurrence  of  an  event 
in  an  incremental  strip  is  independent  of  all  occurrences 
outside  the  strip.   Hence  the  projections  onto  the  coordinate 
axes  give  rise  to  independent  counts  along  the  axes. 

P  (X.Y*)  .  UY*x)VxP(-AY«x)     n  =  Q- 
n  n*  0  <  x  <  X* 

and 

n 


P  (x*,v)  =  ax*yrexp(-xx^  >   n .  0>1>>>< 

0  <  y  <  Y* 


n"~  '"  '  n! 


which  gives  the  univariate  Poisson  distributions  with 
parameters  XY*  and  XX*  respectively. 

Note  here  the  inherent  dependence  on  the  shape  of  the 
assumed  regions.   In  using  rectangular  regions  equal  lengths 
in  the  marginals  reflect  equal  areas  in  the  bivariate 
distribution. 
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If  the  regions  were  circular  then  vertical  projections 
onto  the  axes  would  represent  decreasing  area  as  the  dis- 
tance from  the  origin  increased.   Since  the  occurrence  of 
events  is  assumed  to  be  proportional  to  the  area  projected, 
an  actual  HPPP  would  induce  a  non-homogeneous  process  on 
the  marginals  due  to  the  distortion  in  the  mapping.   For 
clarification,  refer  to  Figure  2.   However,  if  the  regions 
are  circular  then  radial  projections  could  be  utilize:'  so 
that  the  event  occurring  at  (xn,y_)  In  Figure  2  is  repre- 
sented in  the  x-marginal  by  an  event  at  x?.   To  define- 

equal  area  projections  in  this  case  the  transformation 

2 
x  ■--»■  x  =  x'  is  made,  in  which  case  a  unit  increase  in  x1 

defines  the  addition  of  a  unit  amount  of  area  to  the  -  z~- 

For  example,  if  a  unit  area  is  generated  by  a  circle 

radius  r  =  1,  then  the  area  enclosed  in  the  ring  of 

1  <_  r  <_  72  is  the  unit  area,  as  is  the  area  in  the  ring 

f2  _<  r  <_   /T,  etc.   In  general,  Yr\   £  r  <_  yn+1   defines  in 

polar  coordinates  a  ring  with  unit  area. 

Returning  to  the  assumption  of  rectangular  regions, 

three  characteristics  of  the  HPPP  are  now  available  which 

can  be  used  as  the  basis  for  testing  a  sample  for  belonging 

to  the  HPPP  description  of  events  in  a  rectangular  region  P. 

(A)  Independence  of  the  x-ordering  from  the  ordering 
on  the  y-components. 

(B)  Univariate  HPP  (homogeneous  Poisson  process)  in 
the  x-marginal  and,  conditionally  on  n  events  in  R,  a 
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Figure  2.   Vertical  and  radial  projections  of  an 
event  to  form  the  marginal  process.   Shaded  regions 
represent  the  deviations  of  projected  areas  arising 
from  the  rectangular  projection  of  circular  areas. 
Thus,  the  shaded  regions  indicate  the  degree  of 
non-homogeneity  induced  by  the  mapping. 
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uniform  distribution  of  the  distances  to  events. 

(C)   Univariate  HPP  in  the  y-marginal  and,  conditionally 
on  n  events  in  R,  a  uniform  distribution  of  the  distances 
to  events. 

Property  (A)  can  be  tested  against  general  alternatives 
using  a  rank  correlation  procedure  (or  Spearman's  correla- 
tion, see  Pearson  and  Hartley  [1966,  Table  44]).   Properties 
(3)  and  (C)  can  be  tested  by  standard  univariate  methods 
as  in  Cox  and  Lewis  [1966]. 

Note  that  in  the  above  discussion  the  interest  lies  in 
the  nature  of  the  process  rather  than  in  specifically 
describing  the  process.   Thus  the  determination  of  the 
parameter  A  of  the  Poisson  process  is  not  a  current  objec- 
tive and  it  can  be  considered  to  be  a  "nuisance"  parameter. 
Hence  the  conditioning  argument  above  and  the  resulting 
independence  of  the  tests  from  the  value  of  the  parameter 
are  justifiable. 

Now  let  a,  be  the  probability  of  a  Type  I  error  gener 
ated  in  testing  for  randomness,  aR  be  the  corresponding 
probability  in  testing  for  HPPP  in  the  x-marginal,  and  a„ 
likewise  for  the  y-marginal.   Then  the  probability  of  not 
falsely  rejecting  the  HPPP  hypothesis  due  to  the  randomness 
test  is  1  -  a.,  etc.   Hence  the  combined  probability  of  not 
falsely  rejecting  HPPP  is  1  -  prob  {type  I  error}  or 

1  -  P(I)  =  (1  -  cA)(l  -  aB)(l  -  ac). 
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Therefore 

P(I)  =  1  -  (1  -  dA)(l  -  oB)(l  -  ac)  (13) 

is  the  probability  of  falsely  rejecting  a  KPPP  hypothesis. 

If  through  physical  considerations  one  of  the  tests  seems 
more  or  less  significant  than  the  others,  the  analyst  can 
choose  the  weightings  to  so  reflect  the  physical  properties. 
Otherwise  the  values  (and  thus  the  tests)  can  be  weighted 
equally.   This  need  for  the  determination  of  weightings  is 
the  inherent  disadvantage  of  a  multi-level  test. 

The  individual  tests  proposed  above  will  be  briefly 

described.   For  the  rank  correlation  test,  consider  each 

x.  from  (x,y)/.>  which  is  ordered  on  the  y-component.   Also 

consider  the  ordered  realizations  along  the  x  axis,  where 

x.  =  x, . x .   Then 
i    (j  ) 

n  ? 

6  E  (i  -  (j),)2 

1  =  1 _ 

re  «  i  -        ,2     r:         ,  (in) 

s  n(n     -  1)      ' 


where  (j).  is  the  position  of  x.  in  the  x-ordered  sequence, 

is  the  rank  correlation  statistic. 

The  exact  distribution  for  r   can  be  approximated  by 

s 

fitting  a  distribution  to  its  moments  as  discussed  by 

Kendall  and  Stuart  [1951,  p.^77].   The  exact  distribution 

of  r   is  tabulated  in  Biometrika  Tables  For  Statisticians 
s  

[1966,  Table  hk ,  p. 23]  for  observed  values  of  n  between  4 
and  10,  and  the  introduction  to  these  tables  gives  approx- 
imations for  10  <  n  <  20  and  for  n  >  20.   For  10  <  n  <  20, 
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r   can  be  treated  as  a  product-moment  correlation  coefficient 

between  normally  distributed  random  variables.   For  n  >  20, 

r  yn-1     is  assumed  to  be  unit  normal. 
s  ' 

In  testing  the  marginal  distribution  for  HPP,  two 

separate  tests  are  proposed.   First,  the  uniform  conditional 

test  is  used  to  test  against  trends  in  the  data.   As  stated 

in  Cox  and  Lewis  [1966,  p.  153],  "If  the  series  has  been 

observed  for  a  fixed  time  t   {length  X*}  and  n  events  occur 

in  (o,t  ){(o.X*)},  then  the  uniform  conditional  test  is 
o 

based  on  the  variables  U/4 n  -  T./t   {=  X, .  x/X*} (1=1, . . .n) 

(i )    10     (i)        5 

conditionally  on  N,   being  equal  to  n."   The  {brackets}  are 

supplied  to  relate  the  material  in  Cox  and  Lewis  [1966]  to 

this  specific  problem,  and  N   =  n  means  the  number  of 

o 
occurrences  observed  is  n,   Note  that   in  the  conditioning 

of  the  realizations  the  "nuisance"  parameters  XX*  and  AY* 

are  eliminated. 

Secondly,  a  test  based  on  the  ordered  inter-event  spac- 

ings  is  used  to  test  Poisson  against  stationary  event 

processes  which  may  be  non-Poisson.   For  this  test,  Durbin's 

modifications  of  the  uniform  conditional  test  is  used  [Cox 

&  Lewis,  1966,  p.  155].   Referring  to  Figure  3,  Durbin's 

modification  describes  a  transformation  from  the  random 

variable  X  to  the  random  variable  T  and  then  to  the  random 

variable  S. 

Let  3*    -  X*  -  X,  n.   If  the  X,.n.  1  -  1,2,. ..,n, 
n+1         (n)  (i ) 

describe  the  "times  to  events"  in  a  Poisson  process,  then 
the  Tf,  i  =  1,2,..., n+1,  are  independent  exponentially 
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(i)  ~  T(i-1) 


Figure  3-   The  generation  of  the  transformed  variables 
S.  from  the  original  process  X,.x. 
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distributed  random  variables  with  parameter  A.   If  the 

T . ' s  are  then  ordered  and  the  S  's  are  generated  as  shown, 

then  the  S.'s  are  independent  exponential  random  variables, 

where  S.  has  the  expectation  l/( (n+2-i) A) .   Also  the  trans- 

formation  S.  =  (n+2-i)S.  defines  independent  identically 

distributed  exponential  random  variables  with  parameter 

A,  and  therefore  X.  =   E   S.  ,  i  =  l,2,...,n  defines  the 

1    J-l   J 
times  to  events  in  a  Poisson  process  with  parameter  A, 

i    i        n   , 
and  U.  =  ^tt  z  S.   is  the  statistic  upon  which  a  new 

1    X  J-l  j 
uniform  conditional  test  is  based. 

Both  tests  should  be  utilized  as  the  uniform  conditional 
test  is  more  powerful  when  testing  for  trends  while  Durbin's 
modification  is  relatively  more  powerful  in  testing  against 
stationary  event  process  alternatives.   However,  these 
tests  are  not  independent  of  each  other  and  thus  cannot 
be  combined  as  in  (13). 

As  an  alternative  to  the  above  procedure,  the  region 
of  concern  may  be  partitioned  into  several  sub-regions  and 

the  number  of  events  in  each  subregion  used  as  a  basis  for 

2 

X  testing.   This  method  is  discussed  by  Kendall  and  Stuart 

[1951,  pp.  57^-5]  who  mention  the  problem  of  choosing  the 
"right"  partition,  adding  "Whether  a  particular  partition 
has  statistical  interest  depends  on  the  purpose  of  the 
analysis".   Due  to  the  underlying  uniformity  of  the  condi- 
tional distribution,  this  problem  reduces  to  the  selection 
of  the  number  of  regions  which  are  then  used  to  form  equal 
area  sub-regions. 
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Another  alternative  to  the  above  testing  procedure  is 
the  evaluation  of  the  sample  product-moment  correlation 
coefficient  under  the  bivariate  uniform  distribution.   The 
procedure  is  discussed  by  Kowalski  [1972],  but  unfortunately 
his  discussion  does  not  address  the  bivariate  uniform 
distribution.   Kowalski  makes  two  points  very  strongly: 
"Firstly,  the  distribution  of  r  (the  sample  product-moment 
correlation  coefficient  undei  non-normal  assumptions)  may 
differ  from  its  normal-theory  form  and,  secondly,  we  may  be 
in  a  situation  in  which  p  is  a  poor  measure  of  association." 
Hence,  if  the  exact  distribution  for  r  under  the  bivariate 
uniform  distribution  were  known,  then  an  exact  test  for  the 
HPPP  (given  the  occurrence  of  n  events)  could  be  devised. 

Durbin  [1970]  has  also  proposed  distance  methods  for 
testing  bivariate  distributions.   The  process  herein  described 
is  well-suited  to  the  methods  Durbin  uses  since  he  first 
transforms  the  observations  so  that  they  occur  uniformly 
on  the  unit  square.   Hence  the  natural  transformation 
xf  =  x/X*  and  y'  =  y/Y*  avoids  the  problem  of  possible  lack 
of  uniqueness  which  is  the  central  objection  to  the  use  of 
distance  methods.   These  methods  allow  the  analyst  to  adopt 
Durbin' s  bivariate  analog  of  the  Kolmogorov-Smirnov  tests. 
The  advantage  of  this  method  is  the  elimination  of  diffi- 
culties concerning  multi-level  tests  and  partitioning 
tests. 
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The  tests  described  in  this  section  are  very  general  in 
nature,  i.e. 

H  :   the  process  is  KPPP    is  tested  against 
H,  :   the  process  is  not  HPPP. 

Hence  the  alternatives  being  tested  against  are  multitudinous 
If  it  is  desired  to  test  a  realization  as  being  from  a  KPPP 
against  a  specific  form  of  departure  from  HPPP,  better  tests 
may  be  defined  based  on  the  nature  of  the  specific  alterna- 
tive.  For  instance,  one  such  departure  could  be  non-homoge- 
neity, i.e.,  where  X  is  not  considered  to  be  constant  but 
rather  a  function  of  location;  this  subject  is  considered 
in  chapter  III.   Another  departure  might  be  in  the  nature 
of  the  process  itself.   For  example,  events  may  occur 
according  to  some  fixed  plan  in  which  case  the  process  is 
deterministic  and  thus  non-Poisson.   A  process  may  develop 
in  which  the  occurrence  of  an  event  prohibits  the  occurre  :e 
of  another  event  for  some  interval  about  itself,  In  which 
case  events  are  not  independent  of  other  events  and  are 
thus  non-Poisson. 

It  must  be  remembered,  however,  that  tests  against 
specific  alternatives  may  ignore  some  features  that  a  more 
general  test  would  detect  and  thus  each  individual  specific 
test  applies  only  to  the  specific  form  of  departure  being 
considered. 

Moreover,  in  all  reasonable  stationary  alternatives, 
it  does  not  seem  possible  to  derive  the  likelihood  function 
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of  the  observations.   One  thus  cannot  derive  exact  tests. 
For  tests  against  specific  alternatives  based  on  distance 
methods,  see  Holgate  [1972].   Tests  based  on  spectra  are 
discussed  by  Bartlett  [1964]. 

C.   SIMULATING  A  HPPP 

Suppose  one  were  concerned  with  searching  for  submarines 
which  are  assumed  to  be  dispersed  in  such  a  manner  that  the 
locations  at  any  moment  are  generated  by  a  HPPP.   If  one 
search  procedure  is  to  be  selected  from  many  proposed  search 
procedures,  then  a  possible  manner  of  comparing  the  effec- 
tiveness of  the  proposed  procedures  is  to  utilize  each  pro- 
cedure against  several  simulated  dispersions.   In  such  a 
simulation,  the  only  "variable"  which  would  be  of  interest 
would  be  the  procedures,  so  all  variables  such  as  detection 
and  classification  parameters,  facilities  available,  etc., 
would  remain  constant.   Another  problem  which  might  be 
considered  would  be  the  effect  of  the  change  of  such  param- 
eters on  the  search  procedure  selected  (i.e.,  a  sensitivity 
analysis  of  the  procedure  to  assumed  operating 
characteristics  and  facilities). 

By  the  initial  remarks  of  Section  B  above  and  the 
statement  of  equation  (12),  several  methods  of  artificially 
generating  realizations  of  a  KPPP  can  be  determined.   These 
methods  may  then  be  utilized  to  simulate  the  HPPP. 

Assume  that  the  parameter  XX*Y*  is  given.   To  select 
the  number  M  of  events  to  be  observed  in  the  region  with 


29 


area  X*Y*,  generate  a  random  number  U  distributed  uniformly 
on  [0,1].   Set  N  =  n  if 


'y1  (XX^-Y*)1exp(-XXn^)         "  (XX^Y*)1exp(-XX^Y^) 
i         i!  i-1       i! 


2 

The  summations  can  be  evaluated  using  either  x  or  Gamma 

Integral  Tables  [Cox  and  Lewis,  1966,  p. 24]: 


Z    K  .  , =  prob  {y0   >  2u} 

.1       ^     tA2n    MJ 


<»  n-1  -v 

=  /  v 


TH^TTT  dv- 


Next,  consider  a  random  variable  X  distributed  uniformly 

over  (0,X*),  denoted  X  %  U(0,X*),  and  another  independent 

random  variable  Y  %   U(0,Y*).   As  realizations  of  each 

random  variable  are  generated,  number  them  chronologically, 

i.e.  in  order  of  appearance.   Generating  n  (as  determined 

above)  such  realizations  of  each  random  variable  yields  2n 

numbers:  x,  ,...,x  ,y -,,...,  y  . 
1*    '  n*  1'     n 

The  final  problem  remaining  Is  to  select  a  scheme  for 
mating  the  x-  and  y-  realizations  to  form  ordered  pairs 
which  will  constitute  the  realization  of  the  HPPP.   A  few 
such  schemes  are  enumerated: 


1.  The  sequence  <  (x.  jYi)>i =-i  forms  a  HPPP. 

2.  If  the  y.  are  ordered  to  form  <y,.  \>-t==-i»  then  the 
sequence  <(x.sy /iO>  forms  a  HPPP. 
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3.  Similarly,  <^x(±)^1^>    forms  a  HPPP. 

4.  Additionally,  any  random  permutation  of  the  x.  in  2, 
the  y   in  3  or  either  random  variable  in  1  can  be  used 
to  form  a  HPPP.   Thus  <^n+1^±*V ,±* )>    forms  a  HPPP,  etc. 
The  goal  of  the  simulation  and  the  purpose  of  simulating 

the  process  as  a  part  of  the  overall  analysis  must  nov;  be 
considered.   If  during  the  simulation  it  is  desired  to 
generate  independent  realizations  of  the  process,  then  each 
iteration  must  involve  a  selection  of  n,  the  drawing  of  2n 
uniform  variates  and  the  mating  of  the  variates  through  some 
scheme  such  as  those  outlined  in  steps  1-4  above.   On  the 
other  hand,  if  it  is  desired  to  utilize  variance  reduction 
techniques,  then  for  any  drawing  of  2n  random  variates 
several  schemes  could  be  used  for  the  mating  process.   Here 
independence  is  lost  immediately  and  this  loss  must  be 
balanced  by  some  gain  elsewhere  in  the  analysis. 

D.   ESTIMATION  AND  TESTING  FOR  THE  PARAMETER  FROM  A  HOMOGE  EOUS 
PLANAR  POISSON  PROCESS  (HPPP) 

If  the  hypothesis  that  the  process  is  HPPP  with  some 

unknown  value  of  the  parameter  X  is  accepted,  one  might 

like  to  obtain  a  point  estimate  or  confidence  interval 

estimate  for  X,  or  to  test  that  the  process  has  some  given 

parameter  X  .   Note  that  the  parameter  X,  which  was  considered 

to  be  a  nuisance  parameter  in  the  previous  section  where 

the  structural  aspects  of  the  process  per  se  were  tested, 

now  specifies  the  process  completely. 
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Since,  as  was  seen  in  Section  B,  it  is  possible  to 
set  up  the  joint  probability  density  function  of  the 
observations  in  a  HPPP,  point  estimation  of  X  can  be 
based  on  the  method  of  maximum  likelihood.   Note,  however, 
that  each  observation  consists  of  a  single  "look"  at  (or 
realization  of)  the  process  rather  than  n  observations  of 
a  single  random  variable.   Since  it  is  a  stochastic  process 
the  observations  are  not  independent  and  identically  dis- 
tributed.  Hence  the  usual  justifications  for  maximum 
likelihood  procedures  are  not  valid;  see  Brown  [1972]  for 
extensions  of  maximum  likelihood  theory  of  estimation  to 
realizations  of  a  Poisson  process. 

Using  the  results  of  Brown  [1972],  suppose 'that  n  HPPP 
events  are  observed  to  occur  in  a  rectangular  region  of 
area  X*Y*.   From  (11),  for  n  >  0, 

L  =  f((x,y)(1),. . . ,(x,y)(n),n;X}  =  Xne"  (16) 

In  L  ■  nlnX  -  XX*Y*,   (0  <  X  <  *\) 

If  n  ft   0,  this  function  is  -°°  at  X  =  0  and  X  =  »  and  since 

^  ln  L  =  jj.  _  x*Y*   the  slope  of  the  function  decreases 
dX       X 

monotonically  from  <=°  to  -X*Y*.   Thus  ln  L  has  a  unique 

maximum  at  the  point  where  -sj =  0.   Setting  this 

derivative  equal  to  zero  yields  a  unique  maximum  likelihood 
point  estimate  for  X  as 
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where  X  is  unbiased  (since  E(X)  =  —jW   =  X)  and  has  variance 
X/X*Y*.   Note  that  as  the  observed  area  X*Y*  becomes  large, 
the  variance  of  the  estimate  becomes  small;  thus,  by 
Chebyshev's  Inequality  [Lamperti,  1966,  p.  20] 


P{|X  -  X|  >  a}  <  Vag  X  ,   (a  >  0) 


and  as  X*Y*  ->  «>, 

P  {  |  X  -  X  |  >  a }  ■*-  0 

A. 

for  all  positive  a  and  hence  X  converges  to  X  in  probability. 

The  latter  statement  is  equivalent  to  the  assertion  that  X 

is  a  consistent  estimator  for  X.   Also  since  the  variance 

of  X  is  X/X*Y*,  X  has  an  estimated  variance  X/X*Y*  =  n/(X*Y*) 

and  an  estimated  standard  error  of  -/n7x*Y*. 

If  n  =  0,  the  above  method  is  not  applicable.   In  this 

case,  it  might  be  preferable  to  give  a  confidence  interval 

estimate  for  X.   Specifically,  a  one-sided  test  alternative 

is  used  to  generate  a  test  for  the  assumed  value  X   ,,  using 
&  null     ° 

as  an  acceptance  region  only  n  =  0.   Intuitively,  X   ,,  will 

be  small  enough  so  that  X   11X*Y*  <  1  (i.e.,  the  expected 

number  of  observed  events  is  less  than  1).   The  hypothesis 

to  be  tested  is  Hn :  X  =  X   ,,  vs.  H, :  X  >  X  ,,,.   Defining 

0       null      1       null 

a  level  of  significance  a  from  (16)  by 

-X    X*Y* 
prob{N  =  0|X  =  Xnull)  =  1  -  a  =  e   nuir 
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the  hypothesis  HQ  is  accepted  at  the  level  a.   Conversely, 

for  any  given  value  of  a,  A   . ..  may  be  determined  by 

'   null   J  J 


~    XnullX*Y*  "  m  (1  -  ex) 


,       -  In  (1  -  a) 
Anull       X*Y* 


where  the  X   , ,  thus  determined  is  the  largest  value  of  X 
that  the  test  will  accept  at  the  a  level,  given  that  n  =  0. 

Returning  to  the  case  of  n  >  1,  in  order  to  test  that 
the  parameter  of  the  process  has  some  given  value  Xn,  assume 
that  n  events  from  a  HPPP  are  observed  in  a  region  of  area 
X*Y*.   The  hypothesis  to  be  tested  is  II  n :  X  =  Xn  against 
the  two-sided  alternative  H_ :  1  /  L  although  one-sided 
alternatives  can  also  be  considered.   Since  N  is  a  random 
variable  taking  on  all  nonnegative  integer  values  with  sor/ 
positive  probability  for  any  Xn,  there  is  always  some 
possibility  of  an  observed  value  of  the  random  variable  N 
(the  observation  being  denoted  n)  falling  outside  any  finice 
range  of  values.   Thus  a  region  (n  ,n  )  must  be  specified 
such  that  if  N  lies  in  the  region  the  hypothesis  HQ  is 
accepted;  otherwise  the  hypothesis  is  rejected.   The  level 
a  of  the  test  is  the  probability,  given  X  =  XQ,  that  N 
falls  outside  the  region  (n  ,n  ). 

Since  the  test  has  been  defined  to  be  two-sided,  the 
level    is  split  into  upper  and  lower  levels  a   and  a 
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so  that  a  =  a   +  a  .   The  procedure  must  consider  values  of 
A  <  AQ  as  well  as  values  of  A  >    AQ.   To  proceed,  it  is 
necessary  to  define 

P+(n+;A0)  =  P{N  >  n+|A  =  AQ}  =  a+  (18) 

(AnX*Y*)Je   ° 
j=n         J 


and 


P_(n  ;AQ)  =  P{N  <  n  |A  =  AQ}  =  a  (19) 

n"   (AnX*Y*)jexp(-X0X*Y*) 
=  i        __u . 

3=0  j! 

Thus,  for  a  given  a  ,  an  n  may  be  determined  such  that  the 
statement  (18)  just  holds.   Also,  for  a  given  a",  a  n~  may 
be  determined  such  that  (19)  just  holds. 

The  null  hypothesis  is  accepted  at  the  a  level  if  the 
observed  value  of  N  falls  between  the  two  prescribed  limits 
(n   >  n~),  where  prob{N  £    (n  ,n  )}  =  a.   Note  that  as 
stated,  the  result  is  indeterminate  since  a,  once  given, 
leads  to  many  values  for  a  and  a   =  a  -  a  which  satisfy 
the  given  a.   The  manner  of  selecting  a  and  a  must  be 
stated.   Arbitrarily  it  may  bedesirable  to  set  a  =  a   =  a/2 
Asymptotically,  this  choice  of  a  symmetric  acceptance  region 
is  reasonable  since  as  n  increases,  the  distribution  of  N 
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is  approaching  the  (symmetric)  normal  distribution.   The 
choice  of  equal  a  and  a"  may  not  be  reasonable,  however, 
for  small  XQX*Y*  since  the  Poisson  distribution  is 
positively  skewed. 

The  statement  prob{N  I    (n~,n+)|X  =  U  =  a  is  the  result 
of  the  test  of  the  hypothesis  HQ :  X   =  A   at  a  given,  fixed 
level  a.   It  is  this  result  from  which  one  must  usually 
draw  conclusions  regarding  specification  of  the  process. 

If  the  information  thus  available,  i.e.  H  is  rejected 
or  accepted  at  the  pre-determined  a  level,  is  deemed 
insufficient  for  the  purposes  of  a  decision  maker  (for 
example)  then  another  possibility  is  that  the  post-analysis 
information  may  be  extended  by  determining  for  each  obser- 
vation the  exact  a,  a  ,  at  which  the  hypothesis  would  have 
been  rejected.   The  decision  maker  is  then  left  with  the 
problem  of  the  determination  of  his  own  level  of  significance, 
possibly  based  on  his  intuitive  grasp  of  the  problem  and 
its  significance  in  a  larger  frame  of  reference.   Once  he 
has  determined  his  preferred  significance  level,  the  hypoth- 
esis is  rejected  or  accepted  at  the  specified  level  by 
comparison  with  a  .   Thus  the  decision  maker  has  gained 
some  influence  over  the  analysis  but  has  had  to  pay  with 
some  time  to  reflect  on  the  problem  at  hand.   Alternatively, 
he  can  use  a   informally  as  a  "goodness  of  fit"  of  the 
hypothesis . 

Using  (18)  and  (19),  the  significance  test  is  defined 
conventionally  [Cox  and  Lewis,  1966,  p.  30]  to  be:  the 
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hypothesis  A  =  AQ  would  be  accepted  at  the  level  of  signi- 
ficance a  in  a  two-sided  equi-tailed  test  if  the  observed 
number  of  events,  n,  is  such  that  n,  when  used  alternatively 
in  (18)  and  (19)  (i.e.,  is  assumed  to  be  one  or  the  other 
of  the  end-points  of  the  acceptance  region) ,  produces  a 
as  a  solution  to 

P{n;A0}  =  2min{P+(n;A0),P_(n;A0)}  =  ae-  (20) 

Note  that  each  observed  value  of  n  generates  a  new  a   for 

b  e 

any  assumed  AQ,  hence  aQ   =  a  (n,  AQ).   For  example, 
P(30;20)  =  .0436,  P(20;20)  =  .7628  and  P(15;20)  =  .1332. 

It  can  be  seen  that  the  fixed  level  procedure  is 
computationally  simpler,  since  for  a  specified  a  and  A  ,  the 
interval  (n  ,n  )  need  only  be  computed  once  while  in  the 
latter  procedure  a  must  be  recomputed  following  each 
observation  of  N. 

The  inverse  of  the  above  approach  which  utilised  the 
two-sided  equi-tailed  test  of  significance  for  a  given 
value  An  leads  to  the  determination  of  confidence  interval 
estimates  of  A.   Given  that  n  events  are  observed,  it  is 
required  to  determine  some  limits  on  the  range  of  A  such  that 
the  true  parameter  value  A*  lies  within  the  stated  limits 
with  a  probability  1  -  a.   That  is,  it  is  required  to 
establish  a  A~(N)  and  a  A  (N)  such  that 

P{X"(N)  <  A*  <  A+(N)|N  =  n}  =  1  -  a.  (21) 
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Using  P{N  <.n|X=X}<_l-a+to  define  a  X+  as  the 
greatest  X  such  that  equality  just  holds  and  similarly 
using  P{N  <  n|X  =  X-}  =  a"  to  define  a  X"  establishes  the 
limits  such  that  (21)  holds.   For  a  proof  of  this,  see 
Brownlee  [1965 3  p.  121].   Note  that  for  each  realization 
of  N,  a  new  ordered  pair  (X~,X  )  is  defined  so  that  the 
ordered  pair  is  a  function  of  a  random  variable  and  hence 
is  itself  a  random  interval.   The  procedure  only  states 
that  for  (1  -  a)  x  100$  of  the  observations  the  true 
parameter  X*  will  lie  within  the  limits  selected.   The 
limits  for  observed  n  from  0  to  50  are  tabulated  [Pearson 
and  Hartley,  1966,  Table  40]. 

For  a  normal  approximation  to  the  confidence  interval, 
Cox  and  Lewis  [1966,  p.  31]  define  the  upper  a  point  of  the 
unit  normal  distribution  as  c  ,  and  give  the  relationship 

prob{-C.   <  ^r-pr    <   c,  >  =  1  -  a, 

|a   (XX*Y*)1/2  -  |a 

the  relationship  being  correct  as  XX*Y*  •*■   ».   The  confidence 
limits  thus  obtained  are,  to  a  second  degree  of  approximation 
using  a  continuity  correction  and  the  estimate  a(X)  =  yn/X*Y*, 

2 

n  +  2*la  ±  cla^ 
2      2 

For  example,  if  50  events  are  observed  from  a  HPPP,  the 
exact  .05  confidence  interval  is  37-11  <_  XX*Y*  <  65.92 
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[Pearson  and  Hartley,  1966,  Table  40]  v/hereas  the  normal 

approximation  gives  37.79  ±   AX*Y*  <_  66.07. 

2 
There  also  exist  x  approximations  to  the  significance 

tests  and  confidence  intervals  [Cox  and  Lewis,  1966,  p.  33; 

Brownlee,  1965,  p.  173]. 
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III.   NON-HOMOGENEOUS  PLANAR  POISSON  PROCESSES  (NHPPP) 

A.   GENERAL  DISCUSSION 

If  the  stochastic  process  described  above  is  generalized 
to  allow  the  probabilistic  structure  of  the  event  process 
to  be  dependent  on  the  location  of  the  events,  a  non- 
homogeneous  planar  process  is  evidenced.   In  the  simplest 
such  case  a  non-homogeneous  planar  Poisson  process  (NHPPP) 
arises  if,  in  the  definition  of  the  Poisson  process  given 
above,  assumption  I  is  modified  to  become 

I".   There  exists  a  positive  finite  function  X(x,y)  >  0. 
Also  note  that  II  is  changed  by  the  fact  that  the  number  of 
events  in  any  region  is  not  only  a  function  of  the  area 
of  the  region,  but  also  depends  on  the  location  of  that 
region  within  the  universe  under  consideration.   Thus  X 
is  now  expressed  as  X  =  X(x,y),  and  assumption  II  becomes 

II".   prob{N(R  )  =  n} 

{A(A.)}n  exp{-A(A.)} 
_ 

where  A(A.)=  .  /  X(x,y)  dxdy  the  symbol  ./   implying  the 
i    Ai  R.± 

integral  over  an  area  and  X(x,y)  is  assumed  to  be  continuous 
over  R.  (with  area  A.)  so  that  the  integral  is  valid. 

Assumption  III  remains  unmodified,  i.e.  events  occur 
independently  of  any  other  event  or  collection  of  events. 

Under  the  additional  assumption  that  X(x,y)  is  continuous 
within  the  region  of  consideration,  the  incremental 
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development  of  Chapter  II  may  be  extended  to  achieve  a 
description  of  the  NHPPP.   Additionally  the  continuity 
assumption  on  X  and  the  definition  of  the  parameter  in  the 
process  as  an  integral  over  X  eliminates  the  difficulties 
of  line  discontinuities,  although  there  may  be  cases  where 
this  is  an  important  component  of  the  problem.   This  problem 
is  not  considered  here. 

Referring  back  to  Figure  1  in  Section  II-A,  consider 
specifically  the  incremental  strip  defining  region  R, .   If 
the  strip  is  divided  into  n  sub-regions  of  equal  area  by 
taking  n  equal  increments  along  the  x  direction  each  of 
length  6x,  then,  under  the  assumptions  on  the  behavior  or 
^(x,y),  the  process  in  the  ith  sub-region  can  be  approximated 
by  a  HPPP  with  parameter  X~=  X(x,y).  where  (x,y).  is  an 
arbitrary  point  in  the  ith  sub-region.   Specifically  (and 
arbitrarily)   the  lower  left  point  is  chosen  for  the 
succeeding  discussion;  thus  the  parameter  for  the  first 
sub-region  has  parameter  X  =  X(0,Y*).   Continuing,  the 
probability  statements  for  occurrence  of  events  become 

PnCXjy)  =  X(0,Y*)6xAy  +  o(5x,Ay)    0  <  X  <  5x, 
P-^Xjy)  =  X(6x,Y*)6xAy  +  o(6x,Ay)   6x  <  x  <  26x, 


PnCx.y)  =  X(j5x,Y*)6xAy  +  o(6x,Ay)    j6x  <  x  <  (j+l)6x, 
where  j  =  0,1,..., n-1,    Y*  £  y  <_  Y*+Ay   and   n6x  =  X*. 
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Since  the  probability  of  more  than  one  event  in  R 
is  o(X*Ay)  the  probability  statements  above  are  additive 
and 

n-1 
prob  {one  event  in  R  }  =   I  X( j 6x,Y* ) 6xAy  +  o(X*Ay). 

In  the  limit  as  n  ■+  »,  by  the  definition  of  an  integral 

X* 
prob  {one  event  in  R  }  =  {  /   A(x,Y*)dx}Ay  +  o(X*Ay). 

1     0 

By  similar  argument, 

Y* 
prob  {one  event  in  RQ}  =  {  /   A(X*,y)dy}  Ax  +  o(Y*Ax) 


2     0 


and 


prob  {one  event  in  R~,}  =  X(X*  ,Y* )  AxAy  +  o(AxAy). 

By  comparison  with  equations  (3),  (*0  and  (5)  the  above 
statements  lead  to  definitions  for  average  parameters  for 
each  of  the  regions  R, ,  R~  and  R~  as 

X* 
X.(X«,Y*)  =  y*      f      Mx,Y*)dx, 

0 

1    Y* 
X.(X*,Y*)  =  yV  /   X(X*,y)dy,  (22) 

d  0 


and 


X3(X*,Y*)  =  A(X*,Y*). 
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Using  these  average  parameters,  the  equations  (3).,  (*0 
and  (5)  are  generalized,  resulting  in  the  following 
statements : 

prob  {one  event  in  R  }  =  X  (X* ,Y* )X*Ay  +  o(X*Ay), 

prob  {one  event  in  R2>  =  A"2 (X* ,Y* )Y*Ax  +  o(Y*Ax),  (23) 

and 

prob  {one  event  in  R~}  ■  T^CX* ,Y*)AxAy  +  o(AxAy). 

Using  the  result  (22)  as  defining  the  parameter  in  each 
of  the  incremental  areas  in  Figure  1,  equations  (6),  (7) 
and  (8")  become 


P  (X*+Ax,Y*)  =  P  (X*,Y*)[l~A0Y*Ax]+P   . (X* ,Y* ) XnY*Ax 
n      '       n    s       2       n-1    '     2 

+  o(Y*Ax),  (24) 


Pn(X*,Y*+Ay)  =  Pn(X*9Y*)[l-X1X*Ay]+Pn_1(X*,Y*)X1X*Ay 

+  o(X*Ay),  (25) 

and 


P  (X*+Ax,Y*+Ay)  =  P  (X*,Y*+Ay)  +  P  (X*+Ax,Y*)  -  P  (X*,Y*) 
n      J         n  n  n 


-  X0AxAy[P  (X*,Y*)  -  P   , (X*,Y»]   (26) 
3      n  n-l 


+  X..X0X*Y*AxAy[P  (X*,Y*)-2P   -,(X»,Y*) 
12       J      n         n-l 


+  Pn_2(X*,Y*)] 

+  o(Y*Ax)  +  o(X*Ay)  +  o(AxAy). 


M3 


Rearranging  terms,  dividing  by  AxAy  and  taking  the 
double  limit  as  Ax  -*•  0  and  Ay  ->  0  yields 

82P  (X*,Y*) 

8xlf =  "  A3[Pn(X*,Y^)  -  Pn_l(X*,Y*)]        (27) 

+  X-,X0X*Y*[P  (X*,Y*)-2P   _(X»Y*)+P   0(X*,Y*)] 
12      n    *      n-1    '     n-2    ' 


which,  together  with  the  boundary  condition  that  P   is  a 
probability  statement,  gives 


P  (y*   v*)    -  (A(X«,Y*))n  exp{-A(X*,Y*)},  n  =  0,1,... 
n^A  ,Y  ;  "     n!  (28) 


where 


X*   Y* 
A(X*,Y*)  =   /    /   X(u,v)  dudv.  (29) 

0    0 


Thus,  the  number  of  events  occurring  in  a  region 
bounded  by  the  coordinate  axes,  x  =  X*  and  y  =  Y*  has  a 
Poisson  distribution  with  mean  given  by  A(X*,Y*).   The 
mean  can  be  considered  to  reflect  the  cumulative  effect  of 
X(x,y)  in  the  region  of  concern. 

If  n  events  from  a  NHPPP  are  observed  to  occur  in  a 
rectangular  region  defined  as  usual  with  area  X*Y*,  and  the 
events  occur  at  (x,y)/.N,  i  =  l,...,n,  the  labelling  done 
on  the  magnitude  of  the  y-component,  then  the  joint  density 
of  the  events  and  the  probability  that  the  number  of  events 
in  X*Y*  is  n  is  given  by 
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n 

n   X((x,y),,,)exp{-A(X*,Y*)}.  (30) 

j  =  l        U; 

Note  that  (30)  is  a  direct  generalization  of  (16). 

Hence  the  NHPPP  can  be  described  in  a  fashion  similar 
to  the  HPPP,  but  the  expressions  have  acquired  increased 
complexity  due  to  the  necessity  for  the  inclusion  of 
integrals  to  define  the  parameters.   The  degree  of  added 
complexity  is  dependent  upon  the  choice  of  the  specific 
functional  form  for  X(x,y).   The  next  section  develops  the 
expressions  for  one  specific  form. 

B.   A  SPECIAL  CASE 

To  consider  the  location  dependent  type  of  process, 
a  particular  form  for  X(x,y)  is  chosen  as 

X(x,y)  =  exp  {a+  Sx  +  yy  +  5xy}.  (31^ 

Note  that  if  0x  +  yy  +  6xy  changes  very  little  ov;_r  the 
range  of  interest  of  x  and  y,  then 

X(x,y)  =  (1  +  6x  +  yy  +  6xy)exp{a}.  (32) 

Other  relationships  may  be  used;  however,  they  may  cause 
necessary  and  untidy  restrictions  on  the  values  which  the 
constants  a,  3,  y  and  6  may  assume.   In  particular,  X(x,y) 
must  be  greater  than  0  and  the  bivariate  exponential 
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polynomial  (31)  ensures  this  with  no  restrictions  on  the 
range  of  the  parameters. 

Additionally,  algebraic  manipulation  of  the  form  reveals 
that  the  curves  of  X(x,y)  =  c,  c  a  constant,  include  some 
interesting  properties. 

1.  If  <5  =  0,  then  In  X(x,y)  =  c  is  a  family  of  straight 
lines  in  the  plane,  intersecting  the  x-axis  at  an  angle 

6  =  tan   (-3/y).   In  this  case  a  clock-wise  rotation  of 
the  coordinate  axes  through  an  angle  6  would  give  an 
exponential  function  of  y  only. 

2.  If  6  ^  0,  then  In  X(x,y)  =  c  describes  a  system  of 
contour  lines  which  form  a  hyperbolic  paraboloid  with  a 
saddlepoint  at  (-y/<S  ,-8/<5 )  as  is  illustrated  in  Figure  4. 
It  may  be  helpful  to  interpret  the  Figure  in  terms  of  a 
section  of  forest  which  has  been  sampled.   The  line  r 
describes  a  possible  direction  of  steepest  ascent  (DSA) 
which  passes  through  or  near  the  region  being  sampled. 
This  DSA  may  not  be  a  topographic  feature  but  rather  a 
mathematical  expression  for  a  possible  increase  in  density 
of  trees  along  some  line.   Obviously,  there  may  exist  a 
strong  correlation  between  this  mathematical  DSA  and  some 
topographic  features.   Note  here  that  along  the  DSA  maximal 
values  for  X(x,y)  are  found  in  the  sense  that  departing  the 
DSA  at  right  angles  leads  to  decreased  values  for  X(x,y), 
i.e.  decreases  in  the  forest  density. 
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-i 


A  decreasing 


Figure  H.      Contour  lines  for  \(x3y)    =  exp{a+6x+yy+5xy } 
Note  asymptotes  and  region  being  described 


(hatched.)   Here   3/y 


8/6   =  4   and 


all  coefficients  are  positive. 


4- 


A 
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3.   The  exponential  form  can  be  extended  with  little  con- 
ceptual difficulty,  but  possibly  greatly  increased  mathe- 
matical difficulty,  to  describe  a  much  wider  range  of 
possible  circumstances.   For  instance,  it  is  reasonable 

to  assume  that  the  DSA  line  will  bend;  hence  terms  such  as 

2         2 
ex  y  and  £xy   and  higher  order  may  be  included  in  the 

exponent . 

For  the  special  form  of  (3D,  the  cumulative  or  integrated 

intensity  function  A(X* ,Y*)  is  given  by  equation  (29) and 

becomes 

A/(X*>Y*}  =  ^xlla-l\/h   *  Ei^Y+6X*)(ffY*)}-El{(Y+6X*)|} 

(33) 
-  ET{(3+6Y*)J}  +  El(^}  , 

where  Ei  (•)  is  the  exponential  integral  and  Ei(  • )  =  c  +  In  ( • ) 

00   (O1 

+  I    ) .  .  ,  where  c  =  .577216  is  a  constant,  as  defined   i 

i=l1'  x 
Jahnke  and  Emde  [19*15,  p.  2], 

The  likelihood  function  for  the  NHPPP  may  be  develop  ;d 

in  a  manner  similar  to  that  used  in  the  discussion  of  the 

HPPP.   The  discussion  leading  up  to  (16)  is  modified  by  the 

fact  that  the  parameter  is  location  dependent  resulting  in 


=  exp{-A(X*,Y*)}   n   X((x,y)M,),   (n>l)       (3*0 

i=l        K±) 


where  (x,y),.v  is  a  labelling  of  the  coordinates  of  the  n 
point  events.   Thus, 
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n 


In  L  =  -A(X*Y*)  +  i   In  X((x,y),,0 

i  =  l  u; 


follows  directly.   For  the  special  case  of  X(x,y)  given  by 
(3D,  equation  (34)  becomes 

n      n      n 
In  L  =  -A(X*Y*)+B  I   x.+y  I   y.+6  Z  x,y,+na         (35) 

i=l  x   i=l  x   i=l  1  1 

where  A(X*Y*)  is  given  by  (33). 

The  above  joint  density,  or  likelihood,  function  pro- 
vides a  functional  form  which  may  be  manipulated  to  accom- 
plish the  two  principal  concerns  of  the  analysis  of  point 
processes:  hypothesis  testing  and  parameter  estimation.   The 
obvious  null  hypothesis  is  H_ :  B  =  y  =  6  =  0,  in  which  case 
the  nonhomogeneous  Poisson  process  is  being  tested  for  homo- 
geneity since  a  non-zero  a  yields  a  constant  parameter 
X  =  exp(a)  >  0.   Should  the  above  null  hypothesis  be  rejected, 
then  the  analysis  proceeds  to  develop  estimates  for  the 
parameters  (3,  Y  and  <5 .   This  phase  of  the  analysis  may 
proceed  differently  depending  on  how  many  and  which  of  the 
parameters  were  tested  as  being  different  from  zero.   The 
complete,  and  most  complicated,  situation  develops  when  all 
parameters  are  determined  to  be  non-zero.   Testing  of 
parameters  is  the  topic  of  Chapter  IV  while  Chapter  V 
discusses  the  estimation  of  parameters  determined  to  be 
non-zero  as  a  result  of  the  testing  procedure. 
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IV.   TESTING  FOR  NON-ZERO  PARAMETERS 

A.   PRELIMINARIES 

It  is  desired  to  formulate  a  method  for  testing  the 
data  (i.e.,  the  number  of  events  and  their  locations)  in 
order  to  determine  which  of  the  parameters  in  the  model 
given  by  (35),  specifically  a,  B,  y  and  6  ,  are  non-zero. 
Note  that  three  assumptions  are  inherent  at  the  outset: 
first,  that  the  NHPPP  model  is  valid;  second,  that  the 
testing  for  homogeneity  in  Section   II-B  led  to  the  rejec- 
tion of  the  hypothesis  of  homogeneity;  and  third,  that  the 
physical  phenomena  can  be  modelled  by  the  NHPPP  given  by 
(3*0  with  the  parameter  X(x,y)  given  by  (3D. 

Testing  the  Poisson  hypothesis  per  se  when  the  function 
X(x,y)  is  not  known  is  a  compound  problem  which  will  not  be 
considered  here.   It  is  analagous  to  the  compound  problem  in 
regression  analysis  of  testing  both  for  an  unknown  regression 
function  and  for  independent  equal  variance  errors. 

From  the  third  assumption,  the  likelihood  function  for 
the  data  is  given  by 

L  =  exp{-A(X*,Y*)Kn  exp{3SXi  +  yEy±  +  aEx^}  ,   (36) 

where,  for  clarity  in  the  future  development,  £  =  exp{a}. 
Conditioning  on  the  occurrence  of  n  events,  n  >_  1  and 
defining  L{ (x,y) qn , . . . , (x,y) ,% |n;X(x,y) }  =  L(n)  leads  to 
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T        n!  exp{BZx   +  y£y   +  6Ex  y  } 
L(n)  = 


prob{N=n}      Y   X 

(  /   /  exp(Bu  +  yv  +  6uv}  dudv) 
0   0 

(37) 

where  L(n)  is  read  "the  likelihood  function  conditioned  on 
the  occurrence  of  n  events."   Note  that  conditioning  on 
the  number  of  events  observed  has  resulted  in  an  expression 
which  is  independent  of  the  parameter  £(or  a),  i.e.  for 
given  B,  Y>  and  6,  n  is  a  sufficient  statistic  for  a.   This 
is  convenient  because  a  here  is  a  "nuisance"  parameter  since 
the  terms  of  interest  are  those  which  would  indicate  non- 
homogeneity  rather  than  the  establishment  of  the  overall 
rate  of  occurrence.   Thus  by  using  the  conditional  likeli- 
hood a  may  be  eliminated  and  the  testing  can  proceed  for 
non-zero  B,  y  and  6.   In  other  words  the  value  of  a  should 
not  influence  the  testing  for  non-homogeneity  parameters. 

If  n  =  0,  certainly  no  departure  from  homogeneity  couM 
be  evidenced  and  hence  this  case  is  covered  by  HPPP;  see 
II-B  above.   Hence  the  case  of  interest  is  n  >^  1. 

Physically,  the  model  (35)  gives  rise  to  a  parameter 
surface  X(x,y)  which  has  the  properties: 

(a)  B  ^  0;   y^O;   6/0:   lnX  forms  a  hyperbolic 

paraboloid  superimposed  on 
a  tilted  plane,  i.e.  some 
"warping"  of  the  tilted 
plane  is  evidenced. 

(b)  B  ¥■   0;   y  ?   0;   6  =  0:   In  X  forms  a  plane,  tilted 

with  respect  to  the  x-y  plane 

(c)  b  =  Y  ■  0;   6/0:   In  X  forms  a  hyperbolic 

paraboloid. 
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(d)   3=y=6=0:  In  A  forms  a  plane  parallel 

to  the  x-y  plane,  i.e.  a 
HPPP  is  evidenced. 

There  are  a  number  of  possibilities  for  testing: 

(a)   A  test  of 


HQ:   3  =  Y  =  6  =  0 

against 

H-,  :   at  least  one  of  the  parameters  3,  y>  <5  ^  0 

is  a  test  for  non-homogeneity  which  is  more  specific  than 
those  in  Section  II-B  and  is  easily  derived  by  likelihood 
ratio  techniques. 

(b)  The  above  test  is  not  of  great  interest;  generally  the 
specific  non-zero  parameter  is  desired  rather  than  just  that 
at  least  one  of  the  three  is  non-zero.   This  leads  to  the 
question  of  selecting  the  significant  subset,  a  problem 
which  is  difficult  and  as  yet  is  unresolved. 

(c)  The  simpler  problem  is  to  assume  an  ordering,  i.e.  that 
if  M  Y  =  Oj  the  process  is  homogeneous  (6  is  then  assumed 
to  be  0)  and  if  6  or  y   is  non-zero  but  6=0,  then  higher 
order  terms  are  assumed  to  be  zero.   However,  if  the  test 
indicates  non-zero  3   or  y   this  may  be  due  to  an  aliasing 
effect  because  of  a  non-zero  6.   If  further  testing  of 
6=0  against  6^0  reveals  6^0,  then  it  may  well  be  that 
the  true  situation  is  3  =  y  =  0  but  6^0.   The  procedure 

to  be  followed  will  not  discriminate  this  case. 
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The  same  aliasing  effect  occurs  in  testing  of  6  =  0 

against  6^0  where  8  and  y  are  non-zero  and  it  is  desirable 

to  perform  this  test  without  the  effects  of  the  non-zero 

8  and  y.      These  are  thus  nuisance  parameters,  as  was  the 

case  with  a  in  testing  8  and  y.      For  the  present  model  (35), 

one  can  eliminate  these  parameters  because  it  is  seen  from 

the  exponential  form  (36)  that  for  any  6,  (n,  Ex.,  Ey .  ) 

is  a  set  of  sufficient  statistics  for  (a,  8,  y) .   Thus 

6  =  0  is  tested  with  some  function  of  Ex.y.  given  n,  Ex. 

i*  i  °      }   1 

and  Ey  .   This  statistic  has  a  distribution  independent 
of  the  parameters  a,  83  Y* 

The  reason  for  basing  the  conditional  test  on  Ex.y. 
is  that  this  is  (conditionally)  a  sufficient  statistic 
for  6. 

B.   SPECIFIC  TESTS 

Assuming  that  some  ordering  exists  on  the  parameters  as 
discussed  in  possibility  (c)  above,  tests  are  performed 
using  the  sufficient  statistics  (n,  Ex.,  Ey .  ,  Ex.y.)  to 
determine  if  any  non-homogeneity  is  evidenced  by  the  data 
(i.e.,  through  the  statistics).   This  testing  is  more 
specific  in  nature  than  the  testing  encountered  in  Section 
II-B  above  due  to  the  selection  of  a  particular  model. 
The  set  of  sufficient  statistics  arises  from  this  choice  of 
a  specific  model  to  use  as  an  alternative  to  homogeneity. 
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The  testing  will  assume  the  following  sequence: 

(i)   Condition  on  n  and  set  6=0.   Test  H, , . s :  6  =  v  =  0 

0(i)   p    ? 

against  H.^^:   M  0  or  y  /  0.   Note  that  it  would  not  be 
informative  to  test  either  B  or  y  as  a  separate  entity  since 
in  the  formulation  of  the  model  B  and  y  are  unique  only  up 
to  an  angle  of  rotation.   That  is,  testing  of  B  and  y 
jointly  amounts  to  the  detection  of  any  tilt  in  In  A(x,y) 
with  respect  to  the  x-y  plane,  regardless  of  the  direction 
of  the  tilt.   Failure  to  reject  HQ(..v  leads  to  the  assumption 
of  homogeneity  due  to  the  assumed  ordering, 
(ii)   Rejection  of  HQ,.v  leads  to  testing  of 

^O(ii)  :   <5  =  0,  -<»  <  B  <  °°  and   -°°  <  y   <   °° 
against 

HlfiiV   6^0;-co<B<00   and   -°°  <  y    <   °°. 

The  test  thus  specifies  y  and  B  as  nuisance  parameters. 
In  this  test  it  is  necessary  to  first  condition  on  n,  Ex 
and  Zy.  to  eliminate  the  nuisance  parameters. 

In  (i),  conditioning  on  n  and  setting  5=0  leads  to 

(By)n  n!  exp(BExi  +  yly±) 


L(n)  = 


(exp{BX)  -l)n  (exp{yY)  -  l)n 


From  this  it  is  seen  that  the  statistics  (Ex.,Ey.)  are 
(conditionally)  jointly  sufficient  statistics  for  B  and  y 


5^ 


Under  HQ(i), 

Ex±/n  ->  N(X*/2,X*2/12n)   and   Ey  ±/n   ->  M(Y*/2  ,Y*2/12n) 

and  the  statistics  are  independent  (see  Section  II-B).   Hence 
the  expression 

Zx./n  -  X*/2  2     Ey./n  -  Y*/2  2 

(   1  — - )   +  (   1  — ) 

X*/  /l2n  YV  /l2n 

2 

is  asymptotically  Xp»   Rejection  or  acceptance  of  Hn^ . % 

is  based  on  the  adherence  of  the  calculated  value  of  this 

2 
sum  to  the  x  distribution   ,  i.e.  K„  is  accepted  if  this 

sum  has  sufficiently  small  values.  Acceptance  of  H0,.s, 
as  stated  earlier,  leads  to  assumption  of  HPPP;  refer  to 
Chapter  II. 

Following  the  rejection  of  Hn,.>  it  is  necessary  to 
proceed  with  testing  of  Hn,..s.   As  can  be  seen  from  an 
examination  of  (37),  the  complexity  of  the  exact  distribu- 
tion following  another  conditioning  argument  (i.e.  on 
n,  Ex.  and  Ey. )  is  prohibitive.   However,  for  large  sample 
sizes  the  conditional  distribution  can  be  approximated  from 
the  fact  that  Ex./n,  Ey./n  and  Ex.y./n,  conditioned  on  n, 
are  jointly  normally  distributed  for  large  n.   Thus  the 
asymptotic  distribution  of  Ex.y./n,  given  n,  Ex./n  and 
Ey./n,  can  be  found  from  normal  theory  multiple  regression 
results. 
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Under  the  assumption  that  $  ~   y   =   6  *  0 ,  the  trivariate 
normal  distribution  which  arises  is  characterized  by  a  vector 
and  a  matrix.   The  vector  (u_)  of  expected  values  and  the 
variance-covariance  matrix  (E_)  are  given  by 

WW        W 


and 


XtY/2^n 


0  X  Y/24n 
Y2/12n  XY2/24n 
XY2/24n    7X2Y2/lMn 


from  which  p,p  =  C  and  p_ -  =  p  _  =  0.65465. 
In  the  model  given  above, 

H0(ii):   6  =  0;   -"  <  6  <  °°;   -oo  <  Y  < 
is  to  be  tested  against 


Hl(ii):   6  *  0;   -00  <  B  <  °°;   -«  <  y  < 


Since  Ex.y.  is  a  sufficient  statistic  for  6  when  n,  Ix. 
irf  i  •   i 

and  Ey.  are  given,  the  test  can  be  based  on  Ex.y..   Its 
asymptotic  (conditional)  normal  distribution  has  mean 


u   and  standard  deviation  a'   given  by 
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Pxy  "  \y  "  5s  »i3  <EVn  ~   X'2>  + 


a 
+  —^  P2.  (Sy  /n  -  Y/2) 

y 


and 


1 
°xy  =  axy  (1  "  pl3  "  p23  )2  • 

Ex. y./n  -  y ' 

Thus  under  HQr..N,  — — ~ -*■  is  distributed  as  a  unit 

xy 
normal  variate  and  Hq/..n  is  accepted  if  this  statistic  has 

sufficiently  small  values.   Failure  to  reject  Hn,..v  would 

imply  that  the  In  X(x,y). plane  is  tilted  with  respect  to 

the  x-y  plane,  but  no  "warping"  is  evidenced. 

The  above  development  relies  heavily  on  asymptotic 

assumptions.   Small  sample  problems  will  be  much  more 

difficult  to  analyze.   Any  point  in  the  above  procedure  wl  ch 

lead  to  rejection  of  any  hypothesis  would  require  the  ana"  rsis 

to  proceed  with  the  estimation  of  the  non-zero  parameters 

This  is  the  subject  of  the  next  chapter. 


57 


V.   ESTIMATION  OF  PARAMETERS 

It  is  desired  to  formulate  a  method  for  estimating  the 
parameters  a,  3,  y   and  6  of  the  non-homogeneous  planar 
Poisson  model  given  in  IV-A  where  it  has  been  established 
that  a  non-homogeneous  process  is  evidenced  by  the  data. 

Taking  the  logarithm  of  the  conditional  likelihood 
function  (37)  results  in 


In 


L(n)  =  In  n!  +  3Ex.  +  yly .  +  6Ex.y   +  n  In  A ,   (38) 


where  A  =  A(X*,Y*)/£.   Point  estimation  of  (a,  6,  y9    6) 
by  the  method  of  maximum  likelihood  uses  the  conditional 
likelihood  function  (38)  to  develop  the  estimates.   See 
Section  II-D  for  comments  regarding  use  of  maximum  likelihood 
in  this  application.   The  solution  to  the  set  of  simultaneous 
equations 

Y*   X* 
Ex.  -  -r  /    /   u  exp(3u  +  yv  +  5uv}  dudv  =  0 
1   A  0    0 

Y*   X* 
Ey.  -  j     f        f     v  exp(6u  +  yv  +   6uv}  dudv  =  0     (39) 
1    A  0    0 

Y*   X* 
Ex.y.  -  t     f        /  uv  exp(3u  +  yv   +  <5uv}  dudv  =  0, 

1   1       A   Q       Q  ) 

if  obtainable,  provides  the  point  estimates  3,  y   and  6. 
Note  that  this  approach  neglects  the  homogeneous  term 
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during  the  estimation  of  the  parameters  giving  rise  to 
non-homogeneity.   The  neglected  parameter  may  be  estimated 
last . 

v\     /\     s*. 

In  order  for  the  solution  (3,  y>  <5 )  to  equations  (39) 
to  describe  a  relative  maximum  to  In  L | n ,  it  is  necessary 
and  sufficient  that  the  matrix  of  second  partial  derivatives 
(Z)   be  negative  definite,  see  Frisch  [1966,  p.  120].   In 
examining  this  matrix  in  the  case  of  (38),  it  is  helpful 
to  define  S(u,v)  =  exp  {$u  +  yv  +  6uv}.   Then  the  function 

/    v    S(u, v) 
s(u,v)  =   VA? 

has  the  properties: 

(a)  s(u,v)  >  0 

Y*   X* 

(b)  /   /   s(u,v)  dudv  =  1 
0   0 

(c)  s(u,v)  is  continuous  on  [0  <_  u  <_  X*,  0  <_  v  <_  Y*]. 

Hence  s(u,v)  is  a  probability  density  function  [Gnedenko, 
1962,  p.  171]. 

Hence  the  matrix  £  can  be  shown  to  have  diagonal 
elements  such  as 


Y*   X*  -  Y*   X*  2 

a   n  =  -  n[  /   /  u  s(u,v)  dudv  -    (    f        f     us(u,v)  dudv)  ] 
11        0    0  0    0 


=  -  n  Var  U. 
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Continuing,  the  result  is  (where  W  is  defined  to  be  the 
W  =  UV) 

-n  Var  U  -n  Cov  (U,V)     -n  Cov  (U,V/)' 

-n  Cov  (U,V)       -n  Var  V        -n  Cov  (V,W) 
•n  Cov  (U,W)       -n  Cov  (V,W)     -n  Var  W 

and  E_  is  revealed  to  be  a  covariance  matrix.   Note  that  the 
condition  for  a  relative  maximum,  i.e.  E_  negative  definite, 

is  independent  of  the  realizations. 

3 

Now  £_  -   -n  I  where  1_   is  the  usual  variance-covariance 

matrix  for  a  tri-variate  distribution.   But  E_  is  positive 
semi-definite  [Gnedenko,  1962,  p.  212],  hence  -E_  is  negative 
semi-definite.   That  each  of  the  principal  minors  has 
non-zero  determinants  remains  to  be  shown. 

By  the  expressions  given  in  Gnedenko  [1966,  p.  212], 
the  covariance  matrix  £_  can  be  seen  to  be  a  Kankel  matrix 
[Gantmacher,  ls  1959,  p.  338].   Hence  if  the  rows  of  Z_  ar^ 
linearly  independent,  then  the  determinant  of  £_  >  0.   But 

also  Var  U  >  0  since  U  is  a  random  variable  and  Var  U  Var  V  - 

2 
Cov  (U,V)  >  0  since  the  case  of  line  discontinuities  has 

been  excluded  (i.e.,  U  cannot  be  a  linear  function  of  V). 

By  the  same  reasoning,  W  is  linearly  independent  of  U  and 

V.   Hence  all  principal  minors  are  greater  than  zero,  hence 

£_  is  positive  definite,  hence  E_  is  negative  definite. 

S\  >N       A- 

Thus  ($,  Y,  <5 )  provides  a_t  least  a  relative  maximum  to 
In  Lin. 
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Ss.  /S         A 


If  it  were  possible  to  determine  that  (3,  Y,  6)  provides 
a  global  maximum  to  In  L|n  in  the  region  of  interest,  then 
conclusions  as  to  uniqueness  of  the  estimator  could  be 
drawn.   Unfortunately,  global  extrema  are  difficult  to 
establish.   Since  the  method  of  estimation  used  was  maximum 
likelihood,  the  estimates  are  consistent.   Questions  of 
biasedness  are  unresolved. 

In  order  to  solve  the  system  of  equations  (39) »    it  is 
necessary  to  determine  initial  values  for  the  parameters  as 
a  starting  point  for  an  iterative  procedure.   The  partial 
differentiation  of  In L  (35)  with  respect  to  the  parameters 
and  setting  these  partials  equal  to  zero  results,  after 
some  algebraic  manipulation,  in 

n  -A(X,Y)  =  0 

ea  eYY(e(3+5Y)X  -1)    eSX  -  ] 
LX±        6     6  L      B+SY  8    J 

Zyi  +  6A~  X  [ 7+6X g ]  "  °     (40) 


vx   v   +  (il±i-)A  +   ,  e*    „      feYrcgX+YY+6XY  eYY-e6X-n 
LXiyi    l   62  JA+  6(B  +  6Y)(y  +  6X)1  6  Le         "e    e    1J 

-  6XY[eBX+YY+6XY  -  1]  -  6X[eyY  -  1]  -  YY[e6X  -  1]}  =0 


If  it  is  assumed  that  the  sum  $X  +  yY  +  6XY  is  small 
(near  zero)  as  well  as  the  individual  terms  in  the  summation 
being  small,  then  the  exponentials  can  be  approximated  by 
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exp{x)  =  1  +  x,  x  near  zero.   Using  the  first  equation 
in  system  (40)  to  give  the  Value  for  A(X*JY«),  i.e. 
.A(X*,Y*)  =  n,  and  the  linear  approximation  in  the  remaining 
terms  gives  the  abbreviated  system: 


Ex.  +  ?■  n  =  0 


£y,  +  t  n  -  o  CD 


5" 


Lxiyi    +   ^2^n  6(3  +  6Yhy+6X)       L~6~  +   3<5X  Y  +    Y<5XY 


-62XY   +    ByXY]      =      0 


The  solution  to  (4l)  provides  the  initial  estimates  for  the 
parameters.   These  estimates  can  then  be  used  in  (39)  or 
(40)  to  search  for  sequentially  closer  and  closer 
approximations  in  a  mathematical  programming  approach. 

Following  the  determination  of  the  estimates  S,  y 
and  6,  £  can  be  determined  from  the  solution  to  the  first 
equation  in  the  set  (40). 

The  determination  of  confidence  intervals  and  levels 
of  significance  is  not  considered. 
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VI.   CONCLUSIONS 

The  procedures  in  Chapters  III-B,  IV  and  V  are  dependent 
on  the  particular  choice  of  parameter  form;  however,  with 
different  forms  the  concept  of  a  non-homogeneous  planar 
Poisson  process  may  be  used  to  describe  a  wide  variety  of 
"randomly"  occurring  phenomena.   The  choice  of  parameters 
which  may  be  used  is  limited  only  by  assumption  I,  i.e. 
positivity.   One  advantage  of  the  method  discussed  herein 
over  previously  proposed  schemes  is  the  fact  that  the  ' 
specific  form  used  admits  the  possibility  of  a  ridge  or 
line  of  maximum  density  to  be  mathematically  specified 
and  estimated. 

Also  there  is  an  attempt  to  describe  the  underlying 
process  that  caused  the  points  to  appear  where  they  did, 
as  opposed  to  using,  for  instance,  the  arc  within  which 
the  most  events  were  observed  as  the  point  estimate  for 
the  direction  of  maximum  increase. 

Further  efforts  in  this  area  include  a  generalization 
into  four  dimensions  (x,y,z,t)  in  order  that  zoological 
as  well  as  botanical  densities  may  be  studied.   Of  especial 
interest  is  the  estimation  of  densities  of  aquatic  life 
and  how  the  observed  density  fluctuates  with  season  and 
with  changes  in  environment.   The  latter  problem  seems  of 
prime  importance  in  evaluating  the  effects  of  anti-pollution 
programs  on  the  fluid  systems  in  which  plants  and  animals 
exist . 
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Another  problem  which  is  closely  related  to  the  above 
is  that  of  imperfect  sampling  and  how  the  estimates  are 
biased  by  sampling  techniques. 

Chapters  III,  IV  and  V  may  be  redefined  in  terms  of 
data  gathered  within  a  circle  about  some  fixed  point, 
especially  with  consideration  of  the  relative  efficiency 
of  this  data  form  referred  to  by  Matern  [I960]. 
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APPENDIX  A:   THE  BIVARIATE  UNIFORM  DISTRIBUTION 

p 
Given  a  region  R  in  E   of  area  A  and  the  fact  that  the 

probability  of  occurrence  of  an  event  in  any  sub-region  R. 

of  area  A.  within  R  is  simply  A. /A,  a  bivariate  uniform 

distribution  is  described.   For  definiteness  assume  the 

region  R  is  rectangular,  so  A  =  X*Y*.   Nov/ 


Prob  (X  <_   x,Y  <  y)  ■  -22—   =  prob  (X  <  x)  Prob  (Y  <  y) 

X*Y* 


for  0  <_  x  <_  X* ,  0  <_   y  <_  Y* ,  in  which  case  it  is  apparent 
that  the  coordinate  axes  define  independently  chosen 
univariate  random  variables. 

Also,' the  density  function  is  immediately 

f(x,y)  =  1/X*Y*     0  <  x  <  X*,   0  <  y  <  Y*  . 

From  the  density  function  the  joint  density  for  n  independent 
bivariate  uniform  random  variables  is 


f((x,y)1,...,(x,y)n,n)  =  l/(X*Y*)n 


4~  V* 

where  (x,y).  denotes  the  i   pair  of  random  variables 
selected.   Now  n  pairs  of  random  variables,  or  more  simply 
n  points  in  the  plane,  can  only  be  ordered  (without 
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replacement)  in  n!  ways,  independent  of  the  ordering  process 
chosen.   Hence,  the  joint  density  function  for  n  ordered 
bivariate  uniform  random  variables  is 

f((x,y)(1),... ,(x,y)(n),n)  =  n!/(X*Y*)n 

where  (x,y),.«,  is  the  i   point  selected  in  the  ordering 
scheme  utilized. 

As  a  specific  example,  consider  the  n  points  to  be 
labelled  with  respect  to  increasing  magnitude  of  the  y- 
component.   Then 

yk  =  y(k)   k  =  l,..:,n   and   (x,y )  (k)  =  (x^y (k) ) . 

If  the  x-components  are  also  ordered,  then  the  set  of 
points  P  =  ( (x/j\ ,y / • \ ) ;   i,j  =  l,...,n)  defines  n  points, 
of  which  n  are  known  to-be  "occupied,"  that  is,  to  descri  e 
an  event.   For  x,,n,  there  exists  some  j  such  tha'u  y/.\ 
gives  the  y-coordinate  value  for  the  event  which  gave  riL3 
to  x,,v.   Similarly,  for  x,~\  there  are 

now  n-1  j's  remaining,  one  of  which  must  correspond  to  the 
event  giving  rise  to  x,^.   Continuing  to  x,  s,  there  can 
only  be  one  j  left  to  be  associated  with  the  last  x-value . 
Thus  there  are  n!  combinations  of  (x,y)/.x,  i  =  l,...,n 
each  having  density  of  1/(X*Y*)   and  so  the  ordered  bivariate 
uniform  density  is  established  as  that  stated  above. 
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